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■ ABSTRACT 

o 

^ , Using one-dimensional hydrodynamic simulations including interstellar heat- 

ed . ing, cooling, and thermal conduction, we investigate nonlinear evolution of gas 

flow across galactic spiral arms. We model the gas as a non-self-gravitating, un- 
^ , magnetized fluid, and follow its interaction with a stellar spiral potential in a local 

On ! frame comoving with the stellar pattern. Initially uniform gas with density uq in 

the range 0.5 cm~^ < < 10 cm~^ rapidly separates into warm and cold phases 



O 



as a result of thermal instability (TI), and also forms a quasi-steady shock that 
prompts phase transitions. After saturation, the flow follows a recurring cycle: 



OO , warm and cold phases in the interarm region are shocked and immediately cool 

^ ' to become a denser cold medium in the arm; post-shock expansion reduces the 



mean density to the unstable regime in the transition zone and TI subsequently 
mediates evolution back into warm and cold interarm phases. For our standard 
5^ i model with uq = 2 cm~^, the gas resides in the dense arm, thermally-unstable 

transition zone, and interarm region for 14%, 22%, 64% of the arm-to-arm cross- 
ing time. These regions occupy 1%, 16%, and 83% of the arm-to-arm distance, 
respectively. Gas at intermediate temperatures (i.e. neither warm stable nor cold 
states) represents ~ 25-30% of the total mass, similar to the fractions estimated 
from H I observations (larger interarm distances could reduce this mass fraction, 
whereas other physical processes associated with star formation could increase 
it). Despite transient features and multiphase structure, the time-averaged shock 
profiles can be matched to that of a diffusive isothermal medium with tempera- 
ture 1,000 K (which is -C T„arm) and "particle" mean free path of Iq = 100 pc. 
Finally, we quantify numerical conductivity associated with translational motion 
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of phase-separated gas on the grid, and show that convergence of numerical re- 
sults requires the numerical conductivity to be comparable to or smaller than the 
physical conductivity. 

Subject headings: galaxies: ISM — instabilities — ISM: kinematics and dynamics 
— methods: numerical — stars: formation 



1. Introduction 

Spiral arms are the most prominent features in disk galaxies. As the interstellar medium 
(ISM) passes through the moderate gravitational potential well of the stellar spiral arms, it 
is strongly compressed and shocked, producing narrow dust lanes in optical images. Active 
star formation is subsequently triggered in high-density clouds inside the arms, resulting in 
downstream optical arms that contain OB associations and giant H II regions distributed in 
a "beads on a string" fashion (e.g., Baade 1963; Elmegreen & Elmegreen 1983; Elmegreen 
et al. 2006; Shetty et al. 2007). Other arm substructures include filamentary gaseous spurs 
(or feathers) seen in optical extinction, IR emission from dust, and Ha emission from star 
formation (e.g. Elmegreen 1980; Scoville & Rector 2001; Scovillc ct al. 2001; Kcnnicutt 2004; 
Willner et al. 2004; La Vignc ct al. 2006; Gordon 2007), and giant molecular associations 
and atomic superclouds seen in CO and H I radio observations (e.g., Elmegreen & Elmegreen 
1983; Vogel, Kulkarni, & Scoville 1988; Rand & Kulkarni 1990; Knapcn ct al. 1993). The 
locations of these arm substructures downstream from the primary dust lanes indicates that 
the shock compression represents the first step in an evolutionary sequence that begins with 
diffuse ISM gas and ends with star formation (for strongly-bound cores) and dispersal (for 
more weakly self-gravitating structures), although it is uncertain whether spiral arms actually 
enhance star formation rate or just organize it (e.g., Gerola & Sciden 1978; Elmegreen & 
Elmegreen 1986; Sleath & Alexander 1996; Seigar & James 2002). 

Studies of galactic spiral shocks date back to Roberts (1969), who used a semi-analytic 
approach to obtain one-dimensional, stationary shock profiles as functions of the distance 
perpendicular to the shocks (see also Fujimoto 1968; Roberts & Yuan 1970; Shu et al. 1973). 
Woodward (1975) used time-dependent calculations to show that spiral shocks in local models 
indeed develop within one or two crossings of the background arm potential. This and sub- 
sequent work (e.g. Kim & Ostriker 2002) suggests spiral arm shocks in the one- dimensional 
approximation are highly stable for a range of the arm strength. On the other hand, spiral 
shocks have been shown to be intrinsically unstable when the vertical dimension is included 
(Martos & Cox 1998; Gomez & Cox 2002, 2004; Boley & Durisen 2006; Kim & Ostriker 2006). 
Since the arm-to-arm crossing periods are in general incommensurable with the vertical oscil- 
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lation periods, the gas streamlines are not closed, giving rise to shock flapping motions that 
dump a significant amount of random kinetic energy in the gas (Kim, Kim, & Ostriker 2006). 
Under certain (strong compression) conditions, two-dimensional in-plane spiral shocks can 
also become unstable due to strong shear within the arm (Wada & Koda 2004; Dobbs & 
Bonnell 2006). However, these in-plane modes are stabilized by moderate magnetic fields 
(Shetty & Ostriker 2006; Dobbs & Price 2008), and suppressed in fully three-dimensional 
models due to vertical dynamics (Kim & Ostriker 2006). 

Inclusion of gaseous self-gravity tends to enhance the arm response and symmetrize the 
density profile (Lubow et al. 1986), and causes the shock front to move downstream relative 
to the minimum in the potential (Kim & Ostriker 2002). High post-shock density enhances 
the growth of self-gravitating perturbations within spiral arms, although postshock flow ex- 
pansion can hmit this growth (Balbus & Cowie 1985; Balbus 1988). Using two-dimensional 
simulations with both self-gravity and magnetic fields, Kim & Ostriker (2002) demonstrated 
that magneto- Jeans instability (in which magnetic tension forces counterbalance the sta- 
bilizing Coriolis forces) leads to the formation of both arm spurs and GMAs/GMCs with 
realistic properties (see also Lynden-Bell 1966; Elmegreen 1994). Subsequent studies includ- 
ing three-dimensional effects (Kim & Ostriker 2006) and global spiral structure (Shetty & 
Ostriker 2006) have confirmed these findings. 

While recent work has improved our understanding of galactic spiral shocks and their 
larger substructures, these studies have oversimplified the ISM thermodynamics, usually 
adopting an isothermal equation of state. This ignores potential consequences of thermal 
instability (Tl) (Field 1965; see also Meerson 1996 for review), which changes an otherwise 
homogeneous ISM to clumpy, multi-phase gas. (e.g.. Field et al. 1969; Heiles 2001; Wolfire 
et al. 2003). In the classical two- phase picture of the ISM, cold dense clouds are in pressure 
equilibrium with warm intercloud media that surround them (Field et al. 1969). Supernovae 
lead to a hot, diffuse third phase (Cox & Smith 1974; McKee & Ostriker 1977), but because 
massive star formation is spatially correlated and much of the hot gas produced is vented 
away, most of the volume remains relatively unaffected (e.g. Ferriere 1998; de Avillez & 
Breitschwerdt 2004). Since the cold clouds and the warm intercloud gas differ in density and 
temperature by about two orders of magnitudes, their respective responses to spiral shocks 
and downstream expansion flows will be much different from the isothermal case. When 
realistic thermal processes are considered, for instance, warm rarefied gas in the interarm 
region can be converted via shocks to cold dense gas in the arm regions. The reverse phase 
transition can then occur downstream for some fraction of the mass, yielding a quasi-steady 
cyclic exchange. 

Shu et al. (1972) were the first to study the effects of gas coohng and heating on spiral 
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shocks. By considering a mixture of the comoving two stable phases and allowing for phase 
transitions, they calculated steady-state shock profiles for both cold and warm phases. How- 
ever, they employed a pressure- density relation, instead of solving the time-dependent energy 
equation, based on the assumption of instantaneous thermal equilibrium; this precluded the 
possibility of unstable-phase gas in their calculations. 

More recent years have seen a few numerical studies of spiral shocks with explicit heating 
and cooling, but most of these suffer from strong numerical diffusion. Baker & Barker (1974) 
argued that allowance for thermal phase changes produces "accretion fronts/waves" instead 
of spiral shocks, in which the inflowing material radiates its energy away. As they mentioned, 
however, this result could be due to large numerical diffusion; we will indeed show below 
that in a moving medium, numerical conductivity can be large enough to suppress TI. On 
the other hand, Tubbs (1980) and Marochnik et al. (1983) have shown that some models 
develop spiral shocks in which phase transitions from warm interarm gas to cold arm clouds 
occur, although insufficient resolution in their models made the cloud sizes and separations 
significantly overestimated and prevented the transition regions from cold to warm phases 
from being resolved. Very recently, Dobbs & Bonnell (2007) and Dobbs & Price (2008) 
studied the effects of the warm phase on a pre-existing cold phase using particle simulations, 
but they did not allow for phase transitions that arc crucial in spiral shocks with Tl. 

In this paper, we initiate a study of galactic spiral shocks subject to ISM heating and 
cooling, using very high-rcsohition numerical hydrodynamic simulations. We consider one- 
dimensional models that represent slices perpendicular to the arm. We ignore the gaseous 
self-gravity and magnetic fields here, deferring studies of these effects to future work. Our 
primary objectives are to determine overall shock structures under TI, to explore where and 
how the transitions among the cold, warm, and unstable phases occur, and to find statistical 
properties such as temperature distributions, mass fractions, and velocity dispersions. 

The remainder of this paper is organized as follows: in §2 we describe the basic equations 
we solve and present our model parameters and numerical methods. In §3 we test our 
numerical code and quantify the diffusion (due to translational motion over the grid) in 
terms of a numerical conductivity. In §4 we address the evolution of cold and warm phases 
as they traverse spiral shocks, and provide statistical measures to quantify the exchange 
cycle that develops. Finally, we summarize our results and discuss their implications in §5. 
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2. Numerical Methods 

2.1. Basic Equations 

We study galactic gas flows and thermodynamic evolution in response to an external 
stellar spiral potential, which is assumed to be tightly wound with a pitch angle i <^ 1, 
and rotating at a constant pattern speed with respect to an inertial frame. For local 
simulations, it is advantageous to set up a frame corotating with the spiral pattern, centered 
at the position {R,(f)) = {Ro,^pt). The local frame is tilted by an angle i relative to the 
radial direction in such a way that the two orthogonal axes correspond to the directions 
perpendicular (x-axis) and parallel (y-axis) to the local arm segment, respectively (Roberts 
1969). We assume that all physical variables depend only on the x-coordinate, while allowing 
nonzero velocity in the y-direction. Since the independent variable in our local models is a 
projection of the position on a streamline in a large-scale flow onto the x-axis, the distance 
on the X-axis divided by sini represents the distance that the flow has traversed in the 
azimuthal direction along the streamline. Therefore, the temporal interval between one arm 
crossing and the next crossing is the same as it would be for a global model. 

In the local arm frame, the background velocity due to galactic rotation is approximately 
given by 

vo = -Ro(^^o - ^p) sinix + [Ro{Qo - ^p) - qo^ox]y, (1) 

where JIq is the angular velocity of gas at Rq in the inertial frame and = — (d In Q/d In i?) | 
is the local shear rate in the background flow in the absence of the spiral potential (Kim & 
Ostriker 2002, 2006). Assuming that the motions induced by the stellar potential are much 
smaller than Roflo, the basic equations of ideal hydrodynamics expanded in the local frame 
read 

^ + (pvt) = 0, (2) 
dvT 1 

— - + vt • Vvr = — VP - go^^o^^oxY - 20o x v - V$ext, (3) 
ot p 

+ Vr-Ve = 4-PV-vr-p/: + V- (/CVT), (4) 

at 7 — 1 

(see Roberts 1969; Shu et al. 1973; Balbus 1988; Kim & Ostriker 2002; Piontek & Ostriker 
2004), where = Vq + v is the total velocity in the local frame, $cxt is the external stellar 
spiral potential, p£(p, T) is the net cooling function, and /C is the thermal conductivity. 
Other symbols have their usual meanings. We adopt an ideal gas law P = (7 — l)e with 
7 = 5/3. 
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For the stellar spiral potential, we consider a simple sinusoidal shape: 

*ext = *spCOS ( ^ ) , (5) 



analogous to a logarithmic potential of Roberts (1969) and Shu et al. (1973). In equation (5), 
$sp denotes the amplitude of the spiral potential, while — 2t:Rq sin i/m is the arm-to-arm 
separation for an m-armed spiral. We take the size of the simulation domain equal to L^] 
since x varies from —Lx/2 to Lx/2 and $sp < 0, $ext attains its minimum at the center 
{x — Qi). We parametrize the spiral arm strength using 

F^^AM), (6) 



smi 



which measures the maximum force due to the spiral potential relative to the the mean 
axisymmetric gravitational force (Roberts 1969). 

The net cooling function per unit volume is given by pL = n(nA[T] — F), where n = 
p/ {pm-n) is the gas number density and p — 1.27 is the the mean molecular weight per 
particle. For the heating and cooling rates of the atomic ISM, we take the fitting formulae 

F = 2.0 X 10"^^ ergs-\ (7) 

m 



_7 /-1.184xl05\ 2 /^-92\ 

= lO^exp + 1.4 X IQ-^VTexp 

^ V T+ 1000 J \ T J 



cm^, 



F 

suggested by Koyama & Inutsuka (2002) (see also Vazquez-Semadeni et al. 2007). Under 
the adopted net cooling curve, the minimum and maximum pressures for the coexistence 
of the classical warm/cold phases in a static equilibrium are Pmin/^s = 1600 Kcm~^ and 
Pmax/ks = 5000 Kcm~^. The corresponding transition temperatures Tmax = 5012 K and 
^min = 185 K define cold {T < Tmm), warm (T > Tmax), and intermediate-temperature 
phases {T^in < T < Tmax)- To resolve the length scales of TI numerically (e.g., Koyama & 
Inutsuka 2004; Piontek & Ostriker 2004), we include a constant value of thermal conductivity 
/Co = 10^ ergs~^ cm~^K~^.^ The associated Field length below which thermal conduction 
erases temperature perturbations completely is defined by 



Xf = 2n 



/CoT 



(9 In A 
dlnT 



(9) 



(Field 1965). In our models, Xp typically amounts to ~ 0.18 pc. 



^While thermal conductivity is proportional to T^^-^ for neutral hydrogen at kinetic temperatures below 
4.5 X 10* K (Parker 1953), we choose for simplicity a fixed value corresponding to thermal equilibrium at 
T = 1500 K for our standard density no = 2 cm~^. 



- 7- 



2.2. Model Parameters & Numerical Methods 

We consider a simulation box in which the gas is initially homogeneous with density 
uq and pressure Pq when a spiral perturbation is absent. Other than thermal processes 
involving TI, overall dynamics and structures of spiral shocks in our models are completely 
characterized by the arm-to-arm distance and the flow speed Vq^ relative to the perturbing 
stellar potential in the x-direction (as well as sini, qo, F, and ftp/fto)- The azimuthal 
wavenumber m of spiral arms is arbitrary, and the box location and the gaseous angular speed 
relative to the arms can then be specified as Rq = mLx/ {2% sin i) and flQ—flp = 27ivqx/ (ttiLx), 
respectively. To achieve the numerical resolution sufficient to resolve the Field length, we 
consider a small box with = 628 pc. For the relative velocity, we choose Vqx = 13 kms~^ 
corresponding to the rotational velocity of Ro^o = 260 km s^^ with a flat rotation curve 
{qo = 1). The corresponding arm-to-arm crossing time is tcross = L^/vox = 4.7 x 10^ yr, 
which we choose as the fiducial time unit in our presentation. For spiral arm parameters, 
we take pattern speed Qp — Jlo/2, pitch angle sini = 0.1, and strength F — 5% in all the 
models. We note that driven by numerical requirement, spiral arms in our models have a 
small separation and thus a short dynamical time, so that some of our numerical results 
(e.g., mass fractions) that depend on the ratio of cooling time to dynamical time may not 
be applicable to spiral arms with a much larger arm-to-arm crossing time. 

To simulate spiral shocks with varying total gas contents, we consider 12 models that 
have the same initial thermal pressure Po/ks but differ in the initial density no; we adopt the 
solar neighborhood value at Po/ks = 3000 Kcm~^ based on the observational (e.g., Ferriere 
2001; Heiles & Troland 2003; Jenkins & Tripp 2007) and theoretical (e.g., Wolfire et al. 
2003) arguments. Note that although some models are out of thermal equilibrium initially, 
they immediately tend towards equilibrium owing to rapid heating and cooling, with the 
equilibrium pressure depending on no- Table 1 lists the model parameters and simulation 
outcomes. Column (1) labels each run; while the models with the prefix SU rapidly undergo 
TI even with F — 0, the SW and SC models would stay warm or cold throughout were it 
not for the spiral perturbations. Column (2) hsts no- Columns (3)-(8) give the width of, 
and time spent in, the arm, transition, and interarm regions, respectively, in each model. 
The mass and volume fractions of the cold, warm, and intermediate-temperature phases are 
given in columns (9)-(14), respectively. We take model SU2 with tiq = 2 cmT^ as our fiducial 
model; at this density, the total surface density would be ~ 12 pc~^ for vertical scale 
height of ~ 100 pc. 

We integrate the time-dependent partial differential equations (2)-(4) using a modified 
version of the Athena code (Gardiner & Stone 2005). Athena implements a single step, direc- 
tionally unsplit Godunov scheme for compressible hydrodynamics in multi-spatial dimensions 
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and allows a variety of spatial reconstruction methods and approximate Riemann solvers; 
the version we use here employs the piecewise linear method (PLM) with the Roe Riemann 
solver.^ We implement the shearing-periodic boundary condition at the x-boundaries (Haw- 
ley, Gammie, & Balbus 1995). Because of the very short cooling time, energy updates from 
net cooling are made implicitly based on Newton-Raphson iteration, while the conduction 
term is solved fully explicitly. For stable and accurate results, we ensure the time step is kept 
smaller than the CFL condition for thermal conduction as well as the cooling/heating time 
(see Piontek & Ostriker 2004). Our standard models employ N — 16, 384 zones, correspond- 
ing to grid spacing of Ax = 0.04 pc, which satisfies the condition Ax < Af/3 for convergence 
of numerical results (Koyama & Inutsuka 2004); we also ran models with different grid sizes 
in order to study the effects of numerical resolution. 



3. Code Test and Numerical Conductivity 

The Athena code we use has been verified on a wide variety of test problems including 
hydrodynamic shock tubes, advection of a square box in a background shear flow, and one- 
dimensional propagation of sound waves in a rotating, shearing medium. For simulations 
involving cooling/heating and conduction terms, it is crucial to check if the numerical scheme 
employed can resolve the length and times scales of the fastest growing Tl modes. This 
test is of particular importance for the current work since the gas in our models is non- 
static, moving in the x-direction at an average speed of vqx = 13 kms~^, so that numerical 
diffusion associated with advection and zone averaging may reduce the growth rates of TI 
at small scales. In this section, we describe the test results of our numerical code on the 
development of TI in static, rotating, and moving media in the absence of spiral potential 
perturbations, and provide a quantitative measure of the numerical diffusion in terms of 
numerical conductivity. 



3.1. Linear Dispersion Relation 

We test our code by comparing the numerical growth rate of a particular TI mode with 
the corresponding analytic prediction. To this end, we derive a dispersion relation for lo- 
cal, axisymmetric TI in a homogeneous medium that is rotating, shearing, and undergoing 



^Although the piecewise parabolic method (PPM) is known to provide in general less diffusive spatial 

reconstruction than the PLM, our experiments have shown that the PPM often produced negative internal 
energy in the presence of strong radiative cooling inside spiral shocks. 
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uniform translational motion with Vq. We linearize equations (2)-(4) (dropping the exter- 
nal potential), assuming plane- wave disturbances oc e"*~**^^, where n and k are the growth 
rate and wavenumber of the disturbances in the a;-direction, respectively. We follow the 
same steps as in Field (1965), except that we include the non-zero background motions and 
the indirect forces arising from galaxy rotation. The resulting dispersion relation for local 
disturbances is given by 



k \ o , o On q, k ( k \ 9 f k \ 

kT-ko + ^]+KUkT + ^] 



+ n^a(kT + + n{K^ + k'^a?) + a 



0,(10) 



where n = n — ikvox is the Doppler-shifted growth rate, a = {^ksT / jimHY^'^ is the adiabatic 
speed of sound, k = {R~^d{R^Q?') / dR\B^Y^'^ — (4 — 2go)^^^f^o is the local epicyclic frequency, 
and kr-i kp, kjc are the wavenumbers defined by 

7(7 -I) f dC \ 7(7 -I) f dC \ a'p 



0? \dlnT J / " a3 \dlnpj^' 7(7-l)/CT' 

(e.g.. Field 1965). In the limit of go = and Vq^ — 0, equation (10) recovers the dispersion 
relation given by Field (1965) for TI in a rigidly-rotating medium with no translational 
motion. The fact that Vqx occurs in the dispersion relation only through n implies that a 
constant translational motion does not change the growth rates of TI if measured in a frame 
moving with vq^, analogous to Galilean invariance in mechanics. 

For isobaric TI to occur, the term in the square brackets in equation (10) must be 
negative. This of course necessitates kx — kp < 0, the Field criterion for isobaric TI (Field 
1965). Thermal conduction and rotation suppress short and long wavelength perturbations 
against TI, respectively, reducing the unstable range of wavelengths to ki < k < k2, where 
ki and k2 are two positive roots of k'^ + {k)c[kT — kp] + /a^)k'^ -\- ^K^kTk/c/a^ — 0. One 
can show that k2 — > {kic[kp — /ct])^^^ = 2n/Xp for weakly- or non-rotating systems (k — > 0), 
while kf — i> ^K^kr/ {o?[kp — kr]) in the limit of vanishingly small conductivity {kjc — > oo). 
Figure 1 plots as various lines sample growth rates n of TI calculated from equation (10) 
for cases /C = = {dotted line), /C = /Cq with f2o = {solid line), and /C = /Cq with 
rio = 130 kms~^ kpc~^ {dashed line). The density, pressure, and shear parameter are taken 
to be no = 2 cm"^, Po/ks = 3000 cm"^, and qo = 1, corresponding to our fiducial model 
SU2. Stabilization of TI by rotation and conduction are apparent at large and small scales, 
causing the wavelengths of the most unstable disturbances to occur at A ~ 3 pc in between 
the cut-off wavelengths. 
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3.2. Numerical Conductivity 

For our code tests, we consider three kinds of models depending on and flp-. (1) a 
static disk {Qq = Qp = 0); (2) a rotating disk without translational motion {Qq = Qp = 
130 kms~^ kpc^^, yielding Vq^ = from eq. [1]) (3) a rotating disk with translational motion 
{Q.0 = 2Q.p ~ 130 kms~^ kpc~^, so that Vq^ ~ 13 kms~^). The same conductivity K, — ICq 
is adopted for all the models. Other parameters including no, Po, and go are the same as 
in model SU2. In each model, we initialize an eigenmode of TI with wavelength A and vary 
the box size to fit it in, while keeping the grid spacing Ax — 0.04 pc fixed. We monitor 
the evolution of the maximum density and measure its growth rate numerically in the linear 
regime. Figure 1 plots the resulting growth rates from the runs as cross, diamond, and plus 
symbols for the first, second, and third types of models, respectively. Evidently, the numerical 
and analytic results are in good agreement for models without the translational motion, 
confirming the performance of our implementation of the heating/cooling and conduction 
terms. While the translational motion of the gas at a level of — 13 km s""^ does not affect 
large-scale modes much, it significantly reduces the numerical growth rates for small-scale 
modes with A < 2 pc. The discrepancy between the numerical and analytic results appear to 
be due to numerical diffusion that causes the thermal energy to spread out as cooling/heating 
regions are advected with the background flow. 

To quantify the strength of numerical heat diffusion in our code, we solve the linear 
dispersion relation (10) for a given set of parameters to find the effective conductivity /Ccfr 
that yields an analytic growth rate equal to the numerical value obtained from the model 
simulation with the same parameters except for the physical conductivity /C = /Cq- We 
then calculate the numerical conductivity from /C„ = JC^s — In the case of models with 
Vqx = 13 kms~^ and A = 2.56 and 1.28 pc shown in Figure 1, for instance, /C„ = 3.4 x 10^ 
and 1.2 x 10^ ergs~^ cm"-*^ K~^, respectively, about 3.4 and 12 times larger than ICq. 

In order to find the parametric dependences of /C„ on the translational velocity, grid size, 
and perturbation wavelength, we ran a suite of models varying v^x from 2.6 to 26 kms^^. 
Ax from 0.02 to 0.16 pc, and A from 0.32 to 10.2 pc. We also explored the cases with 
different levels of physical conductivity at /C = 0, /Co, or 4/Co. The numerical conductivity is 
calculated only when the numerical growth rate is lower than the analytic value at the same 
/C by more than 1%. Figure 2 plots as various symbols the resulting numerical conductivity 
based on the PLM scheme for spatial reconstruction, showing that /C„ is well fitted by 

This suggests that the numerical conductivity in a moving medium depends rather sensitively 
on the numerical resolution. Equation (12) states that our models with Ax = 0.04 pc 
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have numerical conductivity comparable to the physical conductivity for the fastest growing 
modes, so that the effect of numerical diffusion on the results presented in this paper is not 
significant. 

We repeated the same calculations using the higher-order PPM reconstruction scheme, 
and found the numerical conductivity is linearly proportional to I'oxAx^A"^. Generalizing 

these results, wc rewrite equation (12) as /C„ oc t'oxAa;(A,T/A)^, where p is the order of the 
spatial reconstruction scheme {p = 2 and 3 for PLM and PPM, respectively). This implies 
that the numerical conductivity can be viewed as the diffusion coefficient (oc Vq^Ax) modified 
by the accuracy of the interpolation scheme used in spatial reconstruction (oc [Ax/A]^). 

4. Nonlinear Simulations 

4.1. Standard Model 

We now study the nonlinear evolution of thermally-unstable gas flows under an imposed 
spiral potential. In this subsection, we focus on model SU2 with no — 2 cm~^. We initially 
apply density perturbations created by a Gaussian random field with a power spectrum 
Ipfep oc k"^/^ for 1 < 2'Kk/Lrc < 128 and zero power for 2'Kk/Lx > 128 in the Fourier 
space, corresponding to a one-dimensional Kolmogorov spectrum for sonic disturbances. 
The standard deviation of the density perturbations is fixed to be 1% in physical space. In 
order to suppress rapid motions of the gas caused by an abrupt introduction of the spiral 
potential, we slowly turn it on, reaching the full level F = 5% at t/t 

cross 

Figure 3 shows density distributions of model SU2 at t/tcross = 0, 1, 2.4, and 4.0 together 
with scatter plots of pressure versus density overlaid on the equilibrium cooling curve. The 
gas initially has T = 1500 K and is thus thermally unstable (Fig. 3a). It evolves rapidly to 
cold (n > Tiinax = 8-6 cm^'^ and T < T^^^ = 185 K) and warm {n < n^riin = 1 cm~'^ and T > 
Tmax = 5012 K) phases within 0.5tcross- The TI soon saturates, and random cloud motions 
aided by epicyclic shaking cause cold clouds to collide and merge together, or split sometimes, 
resulting in, on average, 70 cold clouds with a mean cloud separation of ~ 0.8 pc, which 
are in rough pressure equilibrium with the surrounding warm gas at P/ks — 1900 Kcm~^ 
(Fig. 3b). At this time, the spiral potential remains weak and most gas in the unstable 
temperature range corresponds to the boundaries of the cold clouds. As the amplitude of 
the potential grows, the gas is gathered toward the potential minimum, forming a spiral shock 
near a: = 0. The shock reaches maximum strength at around t/tcross = 2.5 shortly after F 
attains the full strength. Both cold and warm phases in the interarm regions continually 
enter the shock front and are compressed to become cold gas with higher density. They 
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subsequently expand and become thermally unstable as they leave the spiral arm regions, 
returning back to the cold and warm interarm phases (see below). At about t/^cross = 3, the 
overall shock structure reaches a quasi-steady state in the sense that the mass and volume 
fractions of each phase do not change appreciably with time, although the shock oscillates 
slightly around an equilibrium position and cold clouds shift as they follow galaxy rotation 
(projected onto the x direction). 

Figure 4 plots the distributions of physical variables in model SU2 at t/tcmss = 4 after 
this quasi-steady state has been reached; many spikes and discontinuities in density as well as 
sawtooth-like velocity profiles arc evident. Figure 5 schematically illustrates the evolutionary 
tracks of the cold and warm phases in the n-P plane. Regions marked with A, B, C, 
and D correspond to interarm, immediate post-shock, spiral arm, and thermally-unstable 
transition zones, respectively; the transition zone refers to the window downstream from the 
arm (between x/L^ ~ 0.03 — 0.2) in Figure 4b, where most of the gas has temperatures 
in the unstable range between Tmin and Tmax- The subscripts 1 and 2 denote the warm 
and cold phases, respectively, of interarm gas, and their respective immediate post-shock 
counterparts. 

In the interarm regions (either x/L^ < or x/L^ > 0.2 in Fig. 4), warm (Ai) and 
cold (A2) phases are in rough equilibrium in terms of the total (= thermal -|- ram) pres- 
sure. While the cold interarm clouds typically have slightly lower thermal pressure than 
the warm interarm medium, their ram pressure is of comparable magnitude, i.e. P/ks ~ 
2000 — 4000 Kcm~^. Some of the cold interarm clouds exhibit "double-horned" structure, a 
consequence of merging with neighbors or splitting into two pieces (see e.g.. Fig. 4a). As the 
warm and cold interarm phases enter the shock front, they experience a strong compression, 
jumping to Bi and B2, respectively. While the density jump is by a factor of 4 for both 
cold and warm phases, the pressure jump for the cold phase is about 10 times larger than 
in the warm phase since the former is about 100 times colder.^ The shocked gas in states 
Bi and B2 is not in thermal balance and subsequently undergoes strong post-shock cooling 
(e.g., Mufson 1974), moving almost isobarically to Ci and C2. The transition from A to C 
is essentially instantaneous (cooling time ~ 10"^ yrs).^ 



•^Strictly speaking, this holds true only for one-dimensional shocks. In two or three dimensions, small 
clouds can experience enhanced compression as the shock wrapping around them is able to propagate towards 
the cloud centers. They may also be subject to dynamical effects such as Kelvin-Helmholtz and Rayleigh- 
Taylor instabilities (e.g., Woodward 1976). 

^Because of its very large post-shock density, the transition from A2 to C2 occurs over an extremely short 
cooling scale 0.02 pc) that is not resolved in our models. 
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One of the characteristics of galactic spiral shocks is that gas accelerates after the max- 
imum shock compression (e.g., Roberts 1969; Balbus 1988; Kim & Ostriker 2002). Because 
of this post-shock expansion, the cold gas inside the arm (0 < x/L^ < 0.03) becomes pro- 
gressively less dense, evolving from C to D. Since the dynamical time scale is much longer 
than the cooling time (due to the reduced velocity inside the arm), the cold gas either at Ci 
or C2 moves all the way down to the transition zone D following the equilibrium curve in 
the n-P plane. When the expanding gas reaches region D (0.03 < x/L^ < 0.2), it becomes 
thermally unstable and turns back into either the warm (Ai) or cold (A2) interarm phase. 
For model SU2, it typically takes about ~ 0.14icross from C to D stages and ~ 0.22icross from 
D to A stages. Gas is in the interarm regime A for the balance of the cycle (~ 0.64icross)- 
When averaged over t/tcross = 5 — 8, the arm, transition, and interarm regions in model SU2 
occupy approximately 1, 16, 83% of the spatial domain, respectively. 

Our models, as well as the real ISM, contain many interarm clouds which would be too 
small to be detected individually in extra-galactic radio observations. In modeling gaseous 
spiral shocks, Lubow et al. (1986) did not explicitly solve for the cloudy structure. Instead, 
they adopted an isothermal equation of state and treated the effects of cold clouds (and 
their collisions) by including a viscous term parameterized by a mean free path Iq for the 
fluid. They found that viscosity renders arm profiles smoother and more symmetric. In 
order to compare our models with a single-phase viscous counterpart, we take a temporal 
average over t/taoss = 5 — 8 of the density distributions in model SU2. We then take a boxcar 
average of the time-averaged profile, with a window of 8 pc. Figure 6 plots the resulting mean 
density profile as a solid line. For comparison, we have run isothermal models that include an 
explicit viscosity term in the momentum equation (3) in a manner similar to in Lubow et al. 
(1986). Selected results are shown in Figure 6. Evidently, larger viscosity tends to produce a 
weaker spiral shock, with a peak that is more symmetric and farther downstream. In terms 
of the strength and placement of the spiral arm, an isothermal model with T = 1000 K and 
Iq = 100 pc provides a fairly good match for the mean density profile when TI is included. 
The isothermal viscous model, however, has a slightly broader arm compared to the time 
average of the multiphase model. Our results thus demonstrate that (for the purposes of 
obtaining a lower-resolution "beam averaged" profile) the effects of clouds embedded in a 
warm medium can be effectively modeled by viscosity (Cowie 1980; Gammie 1996), but only 
if the medium's temperature and the mean free path are chosen appropriately. We note that 
if a temperature comparable to that of the warm medium were adopted for an isothermal 
counterpart, the shock would generally be too weak. 
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4.2. Effects of Initial Number Density 

All the SU models we consider start from a thermally-unstable initial state. Since the TI 
growth time is short compared to the time over which we turn on the spiral potential, these 
models first evolve into a thermally-bistable state before developing shocks. Figure 7 plots 
a density profile at t/tcross = 4 in model SU5 with no — 5 cm~^. Comparing to model SU2 
as shown in Figure 3, Figure 7 shows that the densities (and temperatures) of the cold and 
warm phases after TI saturates are insensitive to the initial gas density, although of course 
models with higher no produce more cold clumps.^ The ensuing development of spiral shocks 
and evolutionary tracks in the n-P plane are also qualitatively similar to those described in 
the previous subsection. 

Perhaps the most notable difference in the late-time states with different no is the size 
of the postshock transition zone. Columns (3)- (8) of Table 1 show that the transition zone 
(and also the arm region, to a lesser extent) widens with increasing Uq. One can easily see 
this by comparing Figures 4d and 7. This trend is of course because models with larger Uq 
have a larger fraction of the total mass in the cold phase that evaporates to maintain an 
unstable state. In addition, the mean separation of cold clouds inside the arm is smaller in 
models with higher no, allowing for merging with neighboring clouds more often during the 
postshock expansion stage. This reduces the expansion rate effectively, thereby extending 
the size of zones with unstable density toward downstream. 

Unlike the SU models, the SW and SC models initially have a low- or high-enough 
density that they are thermally stable and the early development of a spiral shock takes 
place in a single- phase medium. Nevertheless, the shock compression and/or post-shock 
expansion can still produce thermally unstable gas if the initial density is not far from the 
unstable range. For example, the shock in model SW0.5 drives the thermal pressure above 
Pmax, and the gas evolves rapidly via TI into a cold stable phase inside the arm. In model 
SCIO, on the other hand, strong post-shock expansion drives the pressure of the initially 
cold gas below Pmin, after which a fraction of the gas expands to become a warm phase. 
Consequently, the resulting structures at late times in models SW0.5 and SCIO are similar 
to those in the SU models. In models with no further away from the unstable values, 
however, changes in the gas pressure due to shock compression and postshock expansion are 
insufficient to trigger phase transitions. Figure 8 plots equilibrium shock profiles in models 



^Mass conservation requires the number of cold clouds Nc to be given approximately by Nc « no^x/ (Inc), 
where / and ric are the mean size and number density of the clouds. From simulations with differing no, we 
empirically foimd that I ^ 0.8 pc {no/2 cm"'^)'^-^, i.e. the mean cloud size (presumably set by merging and 
splitting), depends only weakly on no- This gives Nc cx ng'^. 
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SWO.l and SC20, which start with uq — 0.1 and 20 cm~^, respectively. The gas in both of 
these models consists only of the warm or cold phase and the profiles are smooth everywhere 
except at the shock front. 

4.3. Temperature Distribution 

Figure 9 plots the volume- weighted and mass- weighted temperature probability distribu- 
tion functions (PDFs) for models SW0.5, SU2, SU5, and SCIO, averaged over t/tcvoss = 5 — 8. 
The vertical dotted lines in each panel indicate T^nm and Tmax, marking the cold, intermediate- 
temperature, and warm phases. For the standard model SU2, the mass-weighed temperature 
PDF is characterized by a broad cold peak at T ~ 20 — 150 K and a narrow warm peak 
at T ~ 6000 — 8000 K, although there also exists a substantial amount of the intermediate- 
temperature gas. In the broad cold peak, the portion of gas with T ^ 50 K corresponds 
to the very dense arm population immediately behind the shock, while the portion with 
50 K ^ T < Tmin includes both arm and interarm cold clouds. Since a larger initial density 
implies a larger initial cold fraction even before the spiral potential is applied, the cold mass 
fractions in both arm and interarm regions increase in models with larger no- For model 
SW0.5, on the other hand, introduction of the spiral potential allows only a small fraction 
of the gas to exceed Pmin and undergo Tl, resulting in a lower cold peak and a higher warm 
peak than in model SU2. For all the models, the cold peak is fractionally much broader (i.e. 
larger AT/T) than the warm peak. 

Of the total intermediate-temperature phase, about 70% is found in the post-shock 
transition zone for models with Tl, while the remaining portion resides in boundary lay- 
ers between the cold and warm phases in the interarm region. Figure 9 shows that the 
distribution of the intermediate phase is almost flat for model SU2 and increasingly fa- 
vors lower temperature as no increases. This is because models with larger no have a 
slower post-shock expansion rate, and thus more gas near Tmin. The mean density and 
density-weighted mean temperature of the intermediate-temperature gas are found to scale 
as 0.41no + 1.29 and Tj 4458/(no + 1.81) K for 0.5 ^ no ^ 10, where rij and no are in 
units of cm~^. For the cold and warm phases, the mean values are ric ~ 23 cm~^, Tc ~ 86 K, 
and nu, ~ 0.39 cm^'^, T^, ^ 6400 K, almost independent of uq. These points lie slightly above 
-Pmin on the cold and warm branches of the thermal equilibrium curve. 
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4.4. Mass and Volume Fractions 

For models with thermal instability, we calculate the mass and volume fractions of 
three phases averaged over t/tcross = 5 — 8. Figure 10 plots against numerical resolution 
the mean mass fractions in model SU2, along with the standard deviations as errorbars. 
While the warm mass fraction is fairly insensitive to the number of grid points N, the cold 
(intermediate-temperature) mass fraction increases (decreases) with increasing resolution for 
N ^ 10^. This is mainly because of the numerical conductivity associated with a large zone 
size for small A^, as explained in §3. When numerical conductivity is large, the boundary 
layers in between the cold and warm phases thicken in proportion to the Field length. This 
increase in the intermediate-temperature mass (in the interarm) occurs primarily at the 
expense of the cold phase (Begelman & McKee 1990; Ferrara & Shchekinov 1993; Inoue et 
al. 2006). In addition, large numerical conductivity tends to suppress TI in the post-shock 
transition zone, yielding less cold gas. 

Overall, the broadening of interarm cloud interfaces in low resolution models is respon- 
sible for 80% of the increase in the intermediate-temperature mass, while the remaining 20% 
is due to the suppression of TI in the post-shock transition zone. There is no interarm cold 
gas in models with A'^ ^ 512, while arm cold gas still exists in these lowest resolution mod- 
els since the post-shock regions are, when fully resolved, denser and broader than interarm 
clouds, and thus are less affected by zone averaging. As long as N ^ 10^, the numerical con- 
ductivity becomes comparable to, or smaller than, the physical conductivity, and the results 
are numerically well resolved. Because most of the difference between moderate and high 
resolution models is at warm/cold interfaces, the extremely high resolutions that we find 
are needed for accurate measurement of the thermally-unstable mass fraction would not be 
required for studies that focus primarily on shock dynamics. We find that the overall shock 
structure (in terms of the breadth of the transition zone, and the time-averaged profiles) are 
comparable to those in the converged models provided N ^ 10^. At these moderate reso- 
lutions, edges of individual clouds are smeared out, but the physically-important transition 
zone downstream from the arm is still recovered. 

In the case of model SU2, the converged values of the mass fractions are 59%, 14%, and 
27% for the cold, warm, and intermediate-temperature phases, respectively. By volume the 
cold takes up 4% of the total, while the warm and intermediate-temperature media occupy 
70% and 26%. Since these proportions are 80%, 10%, and 10% by mass and 10%, 80%, and 
10% by volume when the spiral potential is weak or absent (e.g.. Fig. 3b), this implies that 
spiral shocks and subsequent post-shock expansion zone are important for populating the 
intermediate-temperature portion of the phase plane. 

Columns (9)-(14) of Table 1 give the converged values of the mass and volume fractions 
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of each phase in models with multiple phases. Figure 11 plots these proportions as functions 
of the initial number density no- Interestingly, the fraction of mass in the intermediate- 
temperature phase has a substantial value of /j ~ 0.28, almost independent of Uq.^ The 
volume fraction of the intermediate phase increases with Hq ^ 6 cm^^ and becomes fiat at 
large riQ. Assuming that the cold, warm, and intermediate-temperature phases in each model 
are represented by the characteristic densities Uc, n^, and rii, respectively, mass conservation 
requires the mass fractions in the cold and warm phases to be fc = (1 — n^/ no) — (1 — n^/ nj)/j 
and fyj — riyj/no — {nw/ni)fi, respectively, where ric/riw ^ 1 is assumed (cf. Piontek & 
Ostriker 2004). Dotted hues in Figure 11 plot these theoretical fc and using the empirical 
results fi — 0.28, Uyj — 0.39 cm~^, and rii — 0.41no + 1.29 cm~^. These estimates are overall 
in good agreement with the simulation results. 



4.5. Velocity Dispersions 

Finally, we quantify the level of random gas motions driven in our models due to TI 
and spiral potential perturbations. Spiral arms produce gas streaming motions that are or- 
dered but vary perpendicular to the shock front. Since streaming velocities are much larger 
in amplitude than the true random motions, care is needed in measuring the latter. In 
addition, while the overall shock profiles reach a quasi-steady state after t/tcross = 5, sev- 
eral other effects make it difficult to measure turbulent amplitudes in multiphase models: 
randomly-placed cold clouds move with varying speeds in the interarm region, new struc- 
tures are continuously developing in the transition zone, and the spiral shocks themselves 
undergo small-amplitude oscillations perpendicular to the shock front. To separate out the 
background streaming from the total velocity as cleanly as possible, we first construct ve- 
locity template profiles (vx) and (vy) for each model, where the angle brackets denote a 
time average over t/tcross = 5 — 8. Dashed lines in Figure 4d,e show sample {vx) and {vy) 
profiles in model SU2. We then calculate the density-weighted velocity dispersions using 
af = J p{vi — {v:^Ydx/ J pdx (with i = x or y) in the arm, interarm, and postshock transi- 
tion zones separately. Figure 12 plots o"j(t) for model SU2, while Figure 13 draws the mean 
values (Tj = {crfy^"^ along with the standard deviations AcTj for models with multi-phase 
spiral shocks. 

Figure 12 shows that the velocity dispersions of the cold phase in the arm region ex- 
hibit large-amplitude temporal fiuctuations, with characteristic periods of ~ 0.4tcross and 



^By running models with differing F (not listed in Table 1), we found that fi « 0.040 + 0.045i^(%) for 
3% < F < 7%, but is insensitive to no for the fixed arm parameters. 
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~ 0.2/;cross for and ay, respectively. These variations of ax and ay are due to small- 
amplitude oscillations offsetting the spiral shock front from its mean position. Since stream- 
ing varies strongly within the arm (Fig. 4d,e), the offset of the shock position leads to 
large differences between instantaneous and mean streaming velocities within the arm; these 
contaminate the measured velocity dispersion at locations near the shock front. During the 
quasi-periodic oscillations of the spiral shock front, ax attains its minimum value when the 
shock is maximally displaced upstream, while ay is smallest when the instantaneous shock 
front coincides with its mean position. Therefore, the local minima in the time series of 
ax and ay correspond to the upper limits to the level of random gas motions inside the 
arm. For model SU2, the average values of the local minima in ax{t) and ay{t) are ~ 2 
and 3 kms~^, respectively, which roughly equals ai — Aai. As Figure 13 shows, ai — Aai 
docs not vary much with no, suggesting that the true velocity dispersion within the arm is 
relatively independent of mean density. 

Unlike their behavior in the arm region, the time-averaged velocity profiles in the in- 
terarm and transition zones are relatively smooth, so that the velocity dispersions as cal- 
culated in these regions are relatively free of streaming motions. As Figures 12 and 13 
show, the velocity dispersions in the interarm and transition zones are similar, amount- 
ing to ax ~ 1.3 kms~^, ay ~ 0.8kms~^ for model SU2, and decreasing slowly with in- 
creasing uq. The ratios of the velocity dispersions in the x- to y-directions for all regions 
are consistent with predictions from epicyclic analysis (e.g., Binney & Tremaine 1987), 
o'y/al — ^^/{AQl) = 1 — q/2. Here, the local shear rate q = 2 — (2 — qo)n/no is mod- 
ified from the background value go = —dlnD,Q/dlnR, due to the constraint of potential 
vorticity conservation (Gammie 1996, 2001; Kim & Ostriker 2001, 2002). For a flat rotation 
curve {go = 1), local shear is reversed inside the arm where the local density exceeds 2no; 
this in turn increases (fy/ax ~ (n/2no)^''^ above 1/V2 = 0.7. On the other hand, the whole 
interarm region and most of the unstable region have n < Uq, and thus the ratio ay/ ax is 
smaller than the prediction for a disk without spiral structure. 



5. Summary &c Discussion 
5.1. Summeu-y 

Galactic spiral shocks in disk galaxies play an important role in structural and chemical 
evolution by forming spiral-arm substructures and triggering star formation. Spiral shocks 
inherently involve large variations in the background density, while cooling and heating 
processes that determine the ISM density and temperature depend rather sensitively on 
the local state. The interplay between these processes may significantly alter the shock 
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strengths and structures, compared to those computed under an isothermal approximation. 
In particular, large-scale compressions and expansions across spiral arms may trigger TI, 
thereby regulating transitions among the different ISM phases. In this paper, we have used 
high-resolution numerical simulations to investigate the dynamics and thermodynamics of 
this highly nonlinear process. Our models include heating and cooling terms appropriate for 
atomic gas explicitly in the energy equation, and thus naturally allow for transitions among 
cold, warm, and intermediate-temperature phases. The current investigation employs a one- 
dimensional model in which all the physical quantities vary only in the in-plane direction 
perpendicular to a local segment of a spiral arm. We allow for gas motions parallel to the 
arm, and include galactic differential rotation. The effects of magnetic fields and self-gravity 
are neglected in the present paper. 

Our main results are summarized as follows: 

1. Background flow over the grid, represented by vqx in our models, may result in 
a significant level of numerical diffusion in finite-difference schemes deahng with the cool- 
ing/heating and conduction terms if the numerical resolution is not high enough. We quantify 
the numerical diffusion of our implementation in terms of an effective numerical conductiv- 
ity, and measure it using growing modes of the thermal instability. By comparing growth 
rates with predictions from the linear dispersion relation including conduction, we find that 
the numerical diffusion in the Athena code we use behaves as /C„ oc {voxAx){Ax/\y, where 
Ax is the grid width, p is the order of spatial reconstruction method {p = 2 for PLM 
and 3 for PPM), and A is the spatial wavelength. For typical values of Vq^ = 13 kms~^. 
Ax = 0.04 pc, and A = 3.5 pc in our models, the numerical conductivity with the PLM 
amounts to /C„ ~ 2.3 x 10^ ergs~^ cm~^ K~^, comparable to the physical conductivity 
/Co = 10^ ergs~^ cm~^ adopted in the current work. 

2. Stellar spiral potential perturbations induce shocks that reach a quasi-steady state 
after a few orbits. The resulting flow contains phase transitions provided the mean gas 
density is in the range 0.5 cm~^ <no<10 cm~^, for the parameters considered in the present 
work. Models with no < 0.1 cm~^ or uq > 20 cm~^ yield single-phase spiral shock profiles. 
We divide the flow into three distinct zones based on thermal regime: arm, interarm, and 
transition. The "arm" refers to the highly-comprcsscd postshock region filled with cold gas 
at T < Tmin — 185 K, while in the "interarm" region far from the shock most of the volume 
is occupied by warm gas with T > T^^x = 5012 K. The "transition" zone corresponds to the 
expanding region downstream from the arm where intermediate-temperature gas (Tmin < 
T < Tmax) undergoes TI. Figure 5 summarizes the evolutionary cycle in the density-pressure 
plane: the warm/cold interarm gas (A1/A2) is shocked (B1/B2) and immediately cools to 
become the denser cold arm gas (C1/C2); this subsequently enters the unstable transition 
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zone (D) and evolves back into the warm and cold interarm phases. For our standard model 
SU2 with mean density uq = 2 cm~^, the duration of the arm, transition, and interarm 
stages are approximately 14%, 22%, and 64% of the cycle, respectively, occupying roughly 
1%, 16%, and 83% of the simulation domain. 

3. At late times, instantaneous profiles in models with TI show many density spikes 
representing cold clouds. Time-averaged profiles are, however, relatively smooth, and the 
density peaks representing the arm are more symmetric than in non- diffusive, isothermal 
models. We find that a viscous isothermal model with T = 1,000 K and mean free path 
of /o = 100 pc yields a similar peak density and arm width to the average profile from 
model SU2. This confirms the notion that for purposes where detailed ISM knowledge is not 
needed, multiphase effects on shocks can be approximately treated via a viscosity modeling 
excursions and coUisions of dense clouds, as suggested by Cowie (1980) and Lubow et al. 
(1986). Simulations such as those we have performed are needed in order to cahbrate the 
viscous/isothermal model parameters, however. 

4. For models with multi-phase spiral shocks, intermediate-temperature gas amounts 
to ~ 0.25 — 0.3 of the total by mass, insensitive to no- Of this, about 70% is found in 
the transition zone, while the remaining 30% lies at interfaces between the cold and warm 
media. This suggests that the postshock expanding flows, an inherent feature of galactic 
spiral structure, is important for producing intermediate-temperature gas. The mean density 
and density-weighted mean temperature of the intermediate-temperature phase are found to 
be rii = 0.41^0 + 1.29 cm'^ and Tj = 4458/((no/l cm"''^) + 1.81) K. The fractions of the cold 
and warm phases are 59% and 14% by mass and 4% and 70% by volume for model SU2, 
respectively, and vary with no according to simple expectations based on mass conservation 
with a prescribed density in each phase. 

5. We find that one-dimensional spiral shocks with multi-phase gas produces non- 
negligible random gas motions. At late times, the gas in both interarm and transition zones 
has typical density- weighted velocity dispersions of ax ~ 1.3 kms~^ and ay ~ 0.8 kms~^ in 
the directions perpendicular and parallel to the spiral arms, respectively. This is trans-sonic 
with respect to the cold medium, and subsonic with respect to the warm medium. The cold 
gas in the arms is estimated to have slightly larger values cr-r ~ 2 kms~^ and (T^^ ~ 3 kms"-*^, 
although true turbulence levels may be slightly lower, because it is difficult to fully subtract 
streaming motions in the arm region. 
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5.2. Discussion 

Since the growth rates of pure TI are maximized at the smallest available wavelengths, 
growth of grid-scale noise dominates numerical simulations if it is not suppressed. To prevent 
this numerical problem, it is customary to include thermal conduction that preferentially sta- 
bilizes small-scale perturbations (e.g., Piontek & Ostriker 2004, 2005), and also broadens the 
transition layer between cold and warm phases (e.g., Begelman & McKee 1990). Koyama & 
Inutsuka (2004) studied TI in an initially static medium and showed that numerical results 
converge only if these conductive interfaces are well resolved. They found this requires the 
cell size to be less than one third of the Field length. In modeling systems with large tur- 
bulent (e.g., Gazol et al. 2005) or ordered (this work) velocity flows over the grid, numerical 
conductivity may be considerable if resolution is inadequate. The convergence study pre- 
sented in §4.4 suggests that the convergence criterion of Koyama & Inutsuka (2004) is valid 
in a moving medium, too, with the condition that the grid spacing must be smaller than the 
effective Field length based on the total (physical -|- numerical) conductivity. 

In earlier work, Tubbs (1980) and Marochnik et al. (1983) performed time-dependent 
simulations of galactic spiral shocks that included cooling and heating, but their models 
appear to suffer from large numerical conductivity associated with insufficient resolution. 
While they showed that phase transitions from warm gas to cold clouds do occur in some 
models at the shock locations, they were unable to resolve the post-shock transition zone. In 
particular, the no = 0.5 cm~^ model shown in Figure 7 of Marochnik et al. (1983) exhibits a 
smooth density profile without any indication of TI although much of the region is occupied 
by gas with density and temperature in the unstable range.''' This is presumably because 
large numerical conductivity (due to low resolution) renders the time and/or length scales of 
TI in their simulations longer than the duration and/or width of the post-shock transition 
zone. Indeed, we find that when we run our own models at low resolution (e.g. with 512 
zones or less for model SCIO), the numerical conductivity is large enough that the gas in 
the transition zone does not undergo TI, and instead smoothly converts to interarm warm 
phase. 

In terms of their time-averaged properties, the overall dynamics and flow characteristics 
of spiral shocks with TI remain similar to those of isothermal spiral shocks (e.g., Roberts 
1969; Shu et al. 1973; Woodward 1975; Kim & Ostriker 2002). The sahent features of spiral 



''The heating and cooHng functions adopted by Marochnik et al. (1983) have the critical densities demar- 
cating the warm, intermediate-temperature, and cold phases about 20 times smaller than those in the present 
work. Thus, their model with no = 0.5 cm~^ corresponds to our model SCIO, which develops multi-phase 
structure. 
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shocks with thermal evolution include the facts that they allow phase transitions at the shock 
front (from warm to cold) and in the post-shock transition zone (from cold to warm), and that 
there are cold clumps in the interarm region. Given the very high density and pressure within 
the arm, cold gas downstream from the shock front would transform to molecular clouds if 
self-gravity and chemical reactions for molecule formation were included. By running SPH 
simulations with separate cold and warm "particles" , Dobbs & Bonnell (2007) indeed found 
that the density of the cold component inside arms is sufficient to form molecules (see also 
Dobbs & Price 2008). These authors did not allow for phase transitions between the cold 
and warm components, however, which is an essential aspect of the evolution (cf. Shu et al. 
1972). 

A post-shock transition zone, where the gas is predominantly in the intermediate- 
temperature range, is a necessary feature of any quasi-steady state. About 40% of the 
atomic gas in the local Milky Way is observed to be CNM (Heiles 2001; Heiles & Troland 
2003), which is consistent with our results for mean density no = 1 cm~^, similar to the 
local Milky Way value. ^ Heiles & Troland (2003) also find that about half of the remaining 
atomic gas is consistent with being thermally- stable WNM, while the balance is in the un- 
stable temperature range 500-5000K. These fractions are also comparable to what wc find 
for uq = 1 cm~^. The presence of thermally-unstable gas has been interpreted as owing to 
dynamical effects which act on timescales comparable to the heating and cooling times. Dy- 
namical processes that have been investigated include magnetorotational instability (Piontek 
& Ostriker 2004, 2005, 2007), coUiding flows (Audit & Hennebelle 2005; Heitsch et al. 2005, 
2006; Vazquez-Semadeni et al. 2006), energy injection from OB stars (Vazquez-Semadeni et 
al. 2000; Gazol et al. 2001), and supernova explosions (Avillez & Breitschwerdt 2005; Mac 
Low et al. 2005). Here, we have shown that even without these additional small-scale energy 
sources, a significant amount of intermediate-temperature gas forms as a natural product of 
large-scale spiral shocks with TI, primarily in the expanding region downstream from the 
arm. 

Many numerical studies of Tl have shown that turbulent amplitudes driven by "pure 
Tl" alone are quite small. For example, Piontek & Ostriker (2004) showed that pure Tl in 
an initially static, uniform medium in thermal equilibrium produces only a modest level of 
random gas motions, cr ~ 0.2 — 0.3kms~^. Kritsuk & Norman (2002a,b) similarly found 
that TI of gas starting from thermal equilibrium, or with time-dependent heating, results 
in only subsonic turbulence, a ~ 0.2kms~^ when the cold gas dominates. The velocity 



^The value of no that produces the numerical results consistent with observations may differ in models 
with different since the fractional size of the transition region is likely to depend on the ratio of cooling 
time to arm-to-arm crossing time. 
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dispersions a ~ 1.5 kms~^ of the gas found in the interarm and transition zones in our 
one-dimensional models are significantly larger than those from pure TI. These velocities 
are much lower than those we found in previous (isothermal) simulations (Kim & Ostriker 
2006; Kim, Kim, & Ostriker 2006) in which the vertical shock structure is resolved, however. 
Evidently, shock flapping offers a more effective means of converting galactic rotation to 
turbulence than simple in-plane oscillations. Work is currently underway to see how TI 
interacts with vertically-resolved spiral shocks - and differential buoyancy of cold and warm 
gas - to generate turbulence and structure in the ISM. 
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Table 1. Summary of model parameters and simulation results 
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Note. — Col. (1): Model name; Col. (2): Initial number density in units of cm ^; Cols. (3)-(8): Width and duration of the arm, transition, and 
interarm regions. Cols. (9)-(14): Mass and volume fractions of the cold, intermediate-temperature, and warm phases 
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Fig. 1. — Analytic growth rate of TI versus perturbation wavelength (eq. [10]) for cases 
/C = fio = (dotted), K, = ICq with Q.q = Q [solid], and K, = ICq with = 130 kms^^ kpc~^ 
{dashed), showing that rotation and conduction stabilize perturbations at large and small 
scales, respectively. Symbols indicate the numerical growth rates measured from test sim- 
ulations, all with /C = /Cq: for a static medium {cross), for a rotating medium without 
translational motion {diamond), and for a medium with both rotational and translational 
motion {plus). Note that numerical diffusion associated with translational motion over the 
grid reduces the growth rates at small scales significantly. See text for details. 
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Fig. 2. — Dependence of numerical conductivity /C„ on the translational velocity vqx, the grid 
size Ax, and the perturbation wavelength A of TI, from the test simulations with physical 
conductivity K, — (cross), K — ICq {diamond), and K — 4/Co {plus). Note that Vq^ is 
expressed in units of kms~^, while Ax and A are in parsecs. The solid line, represented by 
equation (12), gives the best fit to the test results based on the second-order PLM for spatial 
reconstruction in the Athena code we use. 
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Fig. 3. — Left: Density snapshots at t/tcross = 0.0, 1.0, 2.4, 4.0 from model SU2. Dotted 
lines labeled by rimax and demarcate the cold and warm phases. Right: Corresponding 
scatter plots of n versus P/kB together with the equilibrium cooling curve {solid line). Dot- 
dashed lines marked by Tmax and T^nin divide the warm (T > T^ax), intermediate-temperature 
{Tmin <T < Tinax), and cold (T < T^i^) phases. 
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Fig. 4. — Profiles of number density, temperature, pressure, and perturbation components 
of the velocities perpendicular {vx — vt,x — vox) and parallel {vy = VT,y — foy) to the arm, 
for model SU2 at t/taoss — 4.0. In (a) to (c), dotted hnes demark the transitions between 
the cold, intermediate-temperature, and warm phases. Dashed lines in (d) and (e) show the 
time-averaged velocity (streaming) profiles over t/tcross = 5 — 8. 
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Fig. 5. — Schematic evolutionary track in the n — P plane {solid lines with arrows) of gas 
cycling through a spiral pattern with TI. The dotted line is the thermal equilibrium curve, 
while dot-dashed lines mark Tmax and Tmin- The warm and cold phases in the interarm 
regions (A) are shocked (B), undergo radiative coohng (C), become thermally unstable due 
to postshock expansion (D), and return to the interarm two-phase state again (A). See text 
for details. 
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Fig. 6. — Density profile averaged over t/tcross = 5 — 8 from model SU2 (solid), compared 
to stationary shock profiles for isothermal models. Dotted curve: T — 1000 K non-viscous 
model; dot-dashed curve: T — 500 K viscous models with Iq — 200 pc; dashed curve: T — 
1000 K viscous model with Iq — 100 pc. 
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Fig. 7. — Density profile of model SU5 with no = 5 cm ^ at t/tcross = 4. Note that the 
postshock transition zone in model SU5 is wider than that of model SU2 shown in Figure 4. 
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Fig. 8. — Density profiles and scatter plots in the density-pressure plane of models SWO.l 
with no — 0.1 cm~^ {upper frames) and SC20 with no — 20 cm~^ {lower frames) at t/tcross — 
2.4. Note that the spiral shocks in these models are smooth and each case contains only a 
single phase of gas. The shock compression is much stronger for all-cold than all-warm gas, 
because of the higher Mach number. 
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Fig. 9. — Mass- weighted {thick lines) and volume- weighted {thin lines) temperature PDFs 
averaged over i/tcross = 5 — 8 for models SWO.l, SU2, SU5, and SCIO. In each panel, vertical 
dotted lines indicate Tmin and Tmax to demark the cold, intermediate-temperature, and warm 
phases. Of the cold component, the gas with T ^ 50 K is located in the arm region, while 
cold clouds with 50 K ^ T ^ Tmin are found in the inter arm or transition zones. 
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Fig. 10. — Mass fractions of the cold, intermediate-temperature, and warm phases in model 
SU2 versus numerical resolution. The numerical results are not affected by numerical con- 
ductivity as long as the number of zones is larger than 10^. 
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Fig. 11. — Mass and volume fractions of the cold (square), warm (triangle), and intermediate- 
temperature (diamond) phases averaged over t/tcmss — 5 — 8, as functions of the initial 
number density Uq. Errorbars indicate the standard deviations of the measurements. Dotted 
lines show the theoretical estimates computed by adopting a constant value /j = 0.28 for the 
intermediate-temperature of the mass fraction. 
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Fig. 12. — Density-weighted velocity dispersions and (Tj^, relative to time-averaged tem- 
plate values, of the gas in the arm (solid), transition (dotted), interarm (dashed) regions of 
model SU2. The large-amphtude fluctuations of the velocity dispersions in the arm region 
are caused by incomplete subtraction of the arm streaming motions. The velocity dispersions 
in both unstable and interarm regions are not subject to this effect, and have mean values 
of ~ 1.3 kms~^ and ~ 0.8 kms~^ in the x- and y-directions, respectively. 
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Fig. 13. — Mean values (symbols) and standard deviations (errorbars) of the density-weighted 
velocity dispersions of each phase during the time span t/tcross = 5 — 8 for models with multi- 
phase spiral shocks. 



